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Abstract 

Huygens' principle states that the solution of the wave equation radi- 
ated by a bounded source can be represented outside the source region 
as a superposition of spherical 'Huygens wavelets' radiated by sec- 
ondary point-sources on a surface enclosing the primary source. This 
was originally proposed as a geometrical explanation of wave propaga- 
tion, but as such it is conceptually problematic because the spherical 
wavelets propagate equally in all directions, thus implying that the 
wave propagates backwards (toward the source) as well as forwards. 
We propose a solution to this problem by generalizing the idea of Huy- 
gens wavelets. Choosing the surface to be the sphere Sr of radius 
R, we show that the Huygens representation of the exterior wave can 
be continued analytically to a complex radius a = R + ia. For any 
real unit vector h, the complex vector an is shown to represent a 
real disk of radius a tangent to 5*^^ at the point Rh. The complex 
sphere Sa consisting of all such vectors ah is therefore equivalent to 
a real tangent disk bundle with base Sr. Just as the points Rn e Sr 
radiate spherical wavelets, so do the tangent disks an G Sa radiate 
well-focused pulscd-bcam wavelets propagating in the outward direc- 
tion n. The analytically continued Huygens formula can be given the 
following real interpretation: the interior wave radiated by the source 
is intercepted by the set of tangent disks ah, which then re-radiate it as 
a set of outgoing pulsed beams. The original wave is thus represented 
in the exterior as a superposition of pulsed beams emanating from 
disks tangent to the sphere Sr, and the coefficients are interpreted as 
local reception amplitudes by the disks. The generalized principle is 
a completeness relation for pulsed-beam wavelets, enabling a pulsed- 
beam representation of radiation fields. Since the new wavelets can be 
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focused by increasing the disk radius a, our construction solves the di- 
rectionality problem of Huygens' original construction. Furthermore, 
it leads to substantial gains in the efficiency of computing radiation 
fields. Only pulsed beams propagating toward the observer need to be 
included and the rest can be ignored while incurring little error. This 
leads to a significantly compressed representation of radiation fields, 
with the compression controlled by the disk radius a. We confirm these 
results by numerical simulations. 
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1 Huygens principle for time-harmonic waves 



Consider a time-harmonic source g of frequency uj supported in a bounded 
volume y C M^: 

q{x, t) = e~'"^^Qu,{x), supp Q^^ C V. 

The wave radiated in free space and observed at the reception event {xr, tr) 
is 

F{xr,tr) = e'''-'^*'''F^{xr) with F^{xr) = j dxG^{xr - x)q^{x), 
where 

G^{r) = , r = \r\ (1) 

r 

is the outgoing fundamental solution of the Helmholtz equation: 

{V^ +uj^)G^{x) = -At:5{x). (2) 

We are using units in which the constant wave propagation speed c = 1, so 
the wave number \s k = lo / c = oj. 

Let be a smooth surface containing V in its interior. Then Green's second 
identity, combined with ([2]) , shows that F^ is given in the exterior of S by 

Fu^iXr) = ^ ^S{X) G^{Xr - x)BMx) (3) 

where dS is the area measure on S, 9„ is the outward normal derivative at 
X G S, and we have introduced the notation 

9{x)<^nfix) = g{x)dnf{x) - dng{x) f{x). 

Equation (3]) is a precise expression of Huygens' principle as formulated by 
Kirchhoff |BC87| IBW99] . It states that in the exterior region, F^{xr) can 
be represented as a superposition of the spherical waves G^{xr — a;), called 
Huygens wavelets, together with their normal derivatives. Hence the points 
X £ S act as secondary sources which collectively form a surface source 
equivalent to the original source g^i in the exterior region]^ 



^The equivalent surface source consists of a single layer {G^jdnf} and a double (dipole) 
layer {-d„Gu: f}. 
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Equation ([3| can be expressed as a condition on the fundamental solution 
by letting be a point source 

with Xe in the interior of S. This gives Fi^{x) = Guj{x — x^), hence ([3| 
becomes 

Gu;{Xr - Xe) = I dS{x)Guj{Xr - x)()nGu;{x - Xe). (4) 

47r Js 

We call Q the Huygens reproducing relation for G^j. To recover (ll), mul- 
tiply by a general source density Quj{xe) supported inside S and integrate 
over Xe- 



Sr 




Figure 1: The sphere Su, the emission and reception points Xe,Xr, and the vectors 



We shall generalize Huygens' principle by continuing analytically in the in- 
tegration variable x, and for this purpose it will be more convenient to work 
with Q than ([s]) . This will be done in the special case where S is the sphere 
of radius R centered at the origin, 

SR = {x = Rh: h£ 5^}, 

where S'^ denotes the unit sphere. Then 

dS{x) = R^dh, where dn = sin^d0d^ 
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is the area measure on 5^, hence Q becomes 

Guj{r) = -— / dhGuj{rr)3j^G^{re), \xe\ < R< \xr\ (5) 

where 



as seen in Figure [T| The normal derivative dn = n - X/ has been replaced by 
the partial derivative dn, and 

G(^(»^e) = , re = re , ^^re = (7) 

re re 

n ( - I I P, - 

'^ijjyfrj — ^ ) ^j. — l^r-l, (JRrr — ^ 

We shall complexify the points Rh of Sn by complexifying and proving 
that this gives an analytic continuation of the distances Ve and r^, hence of 
the right side in ([5]). In the next section we show that this procedure has a 
surprising and beautiful geometric interpretation in real space. 



2 The complex sphere as a tangent disk bundle 

Let a = R + ia ^ C with a > 0, and consider the complexifications of the 
vectors ([6|, 

Ze = an — Xe = Te + lafi, Zr = Xr — ah = Tr — iah, (8) 

regarded as analytic functions of a. To continue Go, {re) and Gi^{rj.) in ([5| 
to C^, we must continue the distances Ve ,rr analytically in a. We will first 
explain the continuation of in detail and then derive the corresponding 
expressions for re- 

The complex distance from ah to Xr is defined by 

Cr = where w = Zr ■ Zr = r'^ — — 2iarr ■ h. (9) 

<^r will be regarded in parallel as an analytic function of 2;,- G and as a 
complex function of Xr G M.^ with ah S fixed. Any analytic function 
f{xr) depending only on Vr can be continued analytically to some domain 
in C'^ by substituting Cr, and we shall regard this as a deformation 

f{Xr) fafii^r) = f{Xr - ah). (10) 
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Since fan is analytic in a, this deformation preserves solutions of differential 
equations such as ([2|J^ The deformation breaks the spherical symmetry of 
Tf. Coupled with a similar deformation of other variables such as r^, this will 
provide a powerful mathematical tool for generating nontrivial and interest- 
ing solutions from simple spherical ones. Furthermore, the singularities of 
deformed solutions give rise to their deformed sources [KSJ. 

Being defined in terms of the complex square root, ("r is double- valued. To 
make it single- valued, a branch cut must be introduced and a branch chosen. 
In the complex variable w G C, we choose the standard branch cut of y/w 
along the negative real axis w <0. But 

w<0 <^ {rr < a, h ■ Vr = 0} <^ {\xr — Rn\ < a, h ■ {xr — Rfi) = 0} , 

hence the branch cut of Cr as a function of Xr G with ah G fixed is 

T>{ah) = {xr : Tr < a, h ■ rr = 0} (11) 
= {Xr : \xr — i?n| < a, h ■ {Xr — Rfi) = 0}. 

This is the disk of radius a centered at Rfi and orthogonal to n, i.e., the 
disk of radius a tangent to the sphere Sfi, at Rn. As a ^ 0, T>(ah) shrinks 
to the one-point set {Rh} and (r — *■ i^'r- We choose the branch with 

Re Cr > 0, so that Cr — > as a — > 0. 

Define the real and imaginary parts of Cr by 

Cr = - illr (12) 

SO that with our choice of branch, 

> and Sgn r]r = Sgn (r^ • h) 

by (|9|. Since w < on V^ah), Cr is imaginary there; hence the branch cut 
can be characterized as 

V{ah) = {Xr : (,r = 0}. (13) 

Choosing cylindrical coordinates {p, (p, z) with the origin at Rh and the z- 
axis along n, ^ and (12) give 

— = £,r — f]1 and ah • = az = ^^ry^ 



^We shall extend this idea to spacetime, where it applies, in particular, to solutions of 
the wave equation and Maxwell's equations. 
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hence 

2 2 2 2 2 2 2/ 2 I a2 2\ a2 2 / 2 , a2\/ 2 2\ 

Thus {S.rj'nr) ^-re related to the cylindrical coordinates (p,z) by 



a/9 = \/a2 + ,f2^a2 - ry^, az = ^rVr- (14) 



This implies the following important inequalities: 

—a < rjr < a and < < '''r (15) 
where the second one follows from Cr — ~ (^^ ~ Vr)- Also by (14), 



+ — r = 1 and ^5 ^ ^ = 1- 



This proves that the level surfaces of and r/j. are 

/>2 7? 



% = ^-^2^ + ^ = U' ^^•>0 (16) 



p2 ^2 



The level surfaces of are the oblate spheroids O^^ and those of rj^ are the 
semi-hyperboloids TCri,.. The restriction zr/^ > follows from az = ^^Vr and 
> 0. As — > 0, Og,. shrinks to the branch disk T){afi) (13). It can be 
shown |K3) that the families O^^ and 7^^^ are mutually orthogonal, forming 
an oblate spheroidal coordinate system deforming the spherical coordinates 
{vr, 6r,(t>r)- They all share a common focal circiej^ which is the boundary of 
the branch disk: 

dV{an) = {xr : r^. = a, n • = 0} = {x.,. : Cr = 0}. (17) 

The last equality shows that d'D{an) is the set of all branch points of (r- 
Whereas f{w) = y/w has a branch point at = 0, C,r{xr) = \/ {xr — an)^ 
has a branch circle. Figure [2] shows &D{afi) and examples of O^^ and TYr/^- 

Since {—Zr)"^ = z"^, is even as a function of Zr G C^. However, it is not 
even as a function of alone. Instead, we have 

(r{—rr — iah) = Cr(''r + ^a"') = Cr(''r — iafi)* . 



^Its physical significance is that it consists entirely of focal points of both O^^ and TLr^^. 
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Figure 2: The real and imaginary parts of (^r — Cr — iVr form an oblate spheroidal 
coordinate system in M.^ centered at Xr = Rfi with the z-axis along h. The third 
coordinate is 0, the standard azimuthal angle. The above plot shows cut-away 
views of an oblate spheroid with = 0.7a, a semi-hyperboloid 'Hr^^ with 
rjr — 0.8a (z > 0) and another with rj^ = — 0.5a(z < 0). Also shown is the focal 
circle d'D{an) with radius a, whose interior is the branch disk 'D{an). 



The last relation is a reality condition or Hermiticity property on the com- 
plex function (^ri^r), and it requires our choice of branch cut ^ 0: 

Cr{z;) = CriZrT- (18) 

We now use this to define the analytic continuation of by 

Ce = V(re + ian)2 = (^\/(re - mn)2^ = (^^ - irj^)* = + irje- 

This shows that Ce and Cr are directed distances. Their sign difference 
indicates that ah is a receiver for the wave propagating from Xe and an 
emitter for the wave propagating to Xr, as illustrated in Figure |3j The sign 
difference is significant because ibD(an) have opposite orientations. 

As functions of x^, for fixed an ^ C^, and rje have the same properties as 
and rjr except that the z-axis is now along — n due to the opposite orientation 
of T>{ah) . For example, the level surfaces of Ce are oblate spheroids and those 
of r]e are semi-hyperboloids with 

Ce > 0, Sgn?7e = Sgn(n • re). 
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Sr 




Figure 3: The complex point ah represents the real disk 'D{an) tangent to the 
sphere Sr at Rn. The complex distances Ce (from to 'D{an)) and (from 'D{ah) 
to Xr) are depicted schematically, emphasizing their directed nature as explained 
in the text. 



However, it is clear from Figure [T] that while r/^ can have any value in [—a, a], 
every emission point X(, in the interior of Sr must have 

n • re > 0, hence < r/e < fl- 
it can be shown that the exact bounds on rj^, as n varies over 5^ are 



It is clear that 7 can depend only on \X(,\ since the minimum of rjf. must be 
spherically symmetric. In particular, 




(19) 



< \Xe\ <R < 7 < 1 



7o = lim 7 



a 



^ < 1, (20) 



and 



Xe = =^ rje = a for all n. 



which is obvious since 



Xe = Ze = an Ce = O =^ = R, "He = d. 
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The functions r]r and r]e will play an important role, and it is helpful to 



interpret them geometrically. From (16) it follows that is asymptotic 



to the cone C^^ making an angle with the positive z-axis and is 
asymptotic to the cone C^^ making an angle i?e with the negative z-axis, 
where 

acosi9r = r]r and a cos -i^e = % • (21) 
See Figure |4} Hence ([l9]) can be restated as 



< i?e < sin~^ ^ (22) 
\a\ 



while < "i^r < TT. 

The parameters ■i?^, 'de are deformations of the spherical coordinates Or, Oe of 
rr,re. A similar interpretation exists for and as deformations of the 
radial coordinates rr,re: the oblate spheroid O^^ containing Xr is tangent to 
the sphere S^^ of radius at the north and south poles, and the same goes 
for O^^ and S^^; that explains why ^ '''r and ^e^fe- These observations 
provide a complete real geometric interpretation of the complex distances 
C,r and Ce in I^^- As an intuitive aid to understanding the idea, think of Cr 
as the 'distance' between the disk D(an) and the point Xy Its complex 
nature reflects the fact that no single real number can characterize this 
distance, and that the distances from Xj. to points on P(an) depend on the 
inclination of the disk, which can be parameterized by 1?^ or ry^. Hence Qr is 
not spherically symmetric, like r^, but cylindrically symmetric around n. 

The functions ^r,ilr simplify if the observer is far from the disk: 



Vr ^ a ^ Cr = \/ r"^ — a'^ — 2iarr ■ n ~ r,. — iafj. • n (23) 
=^ '?r ~ ^r, rjr ^ a COS 6^ where cos Or = Vr ■ n. 

In particular, note that "Or ~ Or as expected. Hence the spheroids O^^ can 
be approximated by the spheres Sr^ and the semi-hyperboloids by their 
asymptotic cones C^^. The deformed variables {^r^'&r) are thus restored to 
their original values (r^. Or). On the other hand, if the observer is far from 
the sphere, then 

\Xr\ ^ R =^ Tr = \Xr — Rn\ ~ \Xr \ — RXr • fl. 

The far-zone approximation assumes that the observer is far from both the 
disk and the sphere, which can be stated succinctly as follows: 

\Xr\ » \a\ =^ Cr = Vl^Crp + — 2aXr ■ h ~ \Xr\ — OiXr ■ fl 
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Figure 4: The semi-hyperboloids Hr]^ (above) and Hjj^ (below) for rjr = 0.9a 
and — 0.5a. Also shown are the asymptotic cones C^^,C^^ and the focal circle 
&D{an). The cones make angles "dr and de with n and —fi, respectively, given by 
(21 1. The Huygens relation in the time domain will favor pulsed beams with rj^ > r]e, 
which means that the wave is focused into narrower propagation hyperboloids Tirj^ 
(asymptotic to the diffraction cones C,?,,) upon being received along 7i^^ and re- 
emitted by 'D{an). 



or 



\xr\ ^ \a\ ~ '^r ~ \^r\ — RcosOr and rjr ~ acosOr- (24) 

In the engineering literature |N86j . the set 

Sa = {aheC^ -.he S^} (25) 
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is called the complex spher^^of radius a in C^. The correspondence 

aheC^ ^ V{ah) C (27) 

establishes a complete equivalence between complex points and real disks 
(where a 'disk' with radius a = is by definition a point). Under this 
equivalence, Sa corresponds to the set of all disks of radius a tangent to Sr, 
which is a tangent disk bundle with base Sr: 

TaiSR) = {V{ah) : n G S^}, a = R + ia. (28) 

3 Generalized principle for time- harmonic waves 

We can now continue ^ to complex space by extending ([7]) to 

Gu,{Ze) = ——, Ze = ah-Xe, Ce = V^e • Ze (29) 
Se 

Guji^Z^^ — — , Zj^ — CC^ O^Tl, (^j, — \/ Zj^ • Zj^. 



If the observer is far from the disk, (23) gives 



> a ^ G^{zr) - 6^'^'=°^^% (30) 

where we have used C,r ^ rr — ia cos Or ~ Vr in the denominator. Thus G^, 
viewed as a function of r,. G M^, has a radiation pattern |HY99] 

For tj > 0, this is the pattern of a beam propagating in the direction of n, 
while for a; < the beam propagates in the direction of —h. The larger waj^ 
the sharper the beam. Note further that these beams are very special in that 



*The term would be more appropriately applied to 

Sc = {z eC^ : z- z = a^}. (26) 

Since Sc, has real dimension 2 for a while Sa has real dimension 4 (complex dimension 
2), Sa is a proper subset of Sa- 

^Recall that c = 1, so Loa = ka — 2iTa/X where k is the wave number and A is the 
wavelength. Thus ua can be interpreted as the number of wavelengths in the circunaference 
of dV^an). 
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they have no sidelobes. That makes them especiahy useful in applications 
such as communications and remote sensing. Analyticity in Zr combines 
with ^ to give 

(V^ + uj^)G^{xr - ah) = when x,. ^ V{ah), (31) 

where is the Laplacian with respect to x^- This proves that the disk 
T>{an) is the source of the beam. Just as the Huygens wavelet Guj{xr — Rh) 
is radiated by a point source 5{xr — Rh) at Xr = Rh, as seen from ([2]), 
so is the beam G^{xr — ah) radiated by the branch disk ^{ah). This 
will be made more precise later, in the time domain. In the limit a ^ 0, 
Guj{xr — ah) becomes the spherical wavelet Guj{xr — Rh). 

This method of deforming spherical time-harmonic waves to beams was first 
introduced by Deschamps |D71] and has become very popular in the engi- 
neering literature under the name complex-source beams, i.e., beams due 
formally to a 'point source' in C^, in our case ah, but interpreted physically 
as a real disk |KS711 IF 761 IC8H IF82] . Solutions to scattering problems where 
the incident field is a complex-source beam are readily obtained by analyti- 
cally continuing solutions with a spherical incident field fCH89j . Complex- 
point receivers were first introduced in [Z SB96] to model directed electroa- 
coustic transducers in ultrasonics, and they have subsequently proven useful 
in cylindrical and spherical near-field scanning for both acoustic and elec- 
tromagnetic fields IHSl EH lH9^ . 

An earlier application of complex distance was made in General Relativity 
by Ted Newman and his collaborators [NJ651 IN65] . who used it to give 
simple derivations of spinning black holes with and without charge (Kerr 
and Kerr-Newman solutions) by deforming known spherically symmetric 
solutions through analytic continuation]^ 

However, none of the above works actually compute the source of G^^. This 
is not trivial because the singularities of G^ on T){ah) are complicated by 
the branch cut: G^ is infinite on the focal circle J- {ah), where Cr = 0, and 
discontinuous on its interior. In |K3] . the source 5^; of G^^ is defined by 



extending p.\ and (31) to 



A-Kd^ixr — ah) = —{y1 + uP')Guj{xr — ah) (32) 



®The first derivation of a cylindrically symmetric solution of Einstein's equation was 
given by Roy Kerr in 1963 |K63| . It was very complicated, which explains why it had taken 
48 years to generalize Karl Schwarzschild's spherical solution. Newman's derivation, based 
on the complex distance, was a model of simplicity. 
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where is the distributional Laplacian with respect to Xr- It is proved 
that 6uj is a generahzed function supported on Xr G T>{ah). In |K4aj . it is 
shown that the analytically continued Coulomb potential 

1 o , 

$ = -, z = x-ia£C'^, C = \/z-z 

generates a reai electromagnetic field H') in the complex-analytic form 

-V$ = E + iH. 

which in turn identifies its source T>{a) as a spinning charged disk whose 
boundary moves at the speed of light. This is the flat-space version of the 
Kerr-Newman black hole, studied from a different viewpoint by Newman 
in |N73j ■ This analysis is generalized to higher dimensions in |K0) . where 
a rigorous connection between solutions of Laplace's equation in R""*"^ and 
the wave equation in M"'^ (Minkowski space with n space dimensions plus 
time) is established, generalizing earlier work by Garabedian |G64] . 

We are now ready to state and prove the analytic Huygens principle for 
time-harmonic waves. 



Theorem 1 For given emission and reception points Xe,Xr with \xe\ < \xr\, 
the Huygens reproducing relation ^ extends analytically to complex R in 
the open set 

A = {a £ C : Re a > \xe\, \a\ < \xr\}, a = R + ia. (33) 
For a £ A, it states that 

f ~ 

G^{xr - Xe) = -— / dhG,^{xr - ah)d^G^{an - x^) (34) 
47r J 



or 



Proof: Write (34 1 as 



dh—^d^—^, r = \xr-Xe\. (35) 
4vr J Qr Ce 



£ = 7^(a), 



where the left side C is independent of a as noted. This reduces to ^ for 
a = R with \xe\ < R < \xr\. The right side TZ{a) is analytic as long as 
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neither Xg nor Xj. belong to any of the branch disks I?(an). But the union 
of all these branch disks is the spherical shell 



5« = IJ 2?(an) = {a; G : R < \x\ < ^/ B? + = |q|}, 



(36) 



Tie 52 



hence x^. must be in the interior of 5^ and x^ in its exterior. This means 
that lZ{a) is analytic in A, and since it is constant on the line segment An M, 
it must be constant throughout A. ■ 



Equation (34) can be interpreted physically as follows: G^{ah — x^) is the 
reception amplitude by the disk T){an) of the wave emitted by Xe, which 
in turn stimulates the emission of the complex-source beam Gaj{xr — otn) 
propagating to x^.. The spherical wave Guj{xr — x^,) from x^, to Xy is thus 
represented as a sum of beams. 



Equation ( 35 ) can be further simplified by letting 



Ce(a) = \/ {an - x^)"^ and Cr(a') = \/(x'r 
with a and a' independent. Then 



^jaj(Cr+Ce) 



(37) 



(38) 



where 



da' a = {da' " Qa] L/, 



Applying the derivatives gives a version of (38) more suitable for numerical 
computations: 

■2 <■ dn 



4vr } CrCe 



«'^(Cr - Ce) - 7^ + ^ 



(39) 



where 



C = 9.Ce = ^^^, C = 9„<. = ^— (40) 



Let us note a symmetry of (38). Since the left side satisfies the Fourier 



transform reality condition /(— u;)* = /(w), so must the right side; thus 



47r 



dn 



(41) 
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The branches defined by Re Ce ^ and Re Cr ^ satisfy the reality 



conditions (18) 



Ce{ar = Ce{a*) and Ua')* = CAa'*), 



hence the right side of (|4T|) is simply ( 38 ) with a ^ a* and a' — > a'* . Since 



the set A (33) is symmetric under complex conjugation, this explains why 



(41) and (38) are consistent. That is, the right side of (38) satisfies the 



extended reality condition 



f{-u;,a*r = fiu;,a). 



(42) 



supported in Sr is 



Equation (34) implies that the field radiated by an arbitrary source Quji^e) 

(43) 



a 



where 



(44) 



is the analytic continuation of the radiated field F^{Rn) from Sr to Sa- 
Fi^{an) can be interpreted as the reception amplitude of the radiation field 
by the disk 'D{an) |ZSB96j . See also Section Is] where this is proved in the 



time domain using a rigorous definition of pulsed-beam sources. Thus (43) 
has a simple physical interpretation: the field radiated by g^^ is intercepted 
by'D{an) and re-radiated by'D{a'n) to give an identical field in the exterior, 
showing that g^j can be replaced by an equivalent source on the tangent disk 



bundle TaiSpi) given in (28). 



Equation ( 43 ) gives the field radiated by g^j as a superposition of the complex- 
source beams G^{xr — ah) with source points ah G Sq.. The first exact 
representation of this type was obtained by Norris |N86] . who expressed 
the field of a single real point source at the origin in terms of complex- 
source beams emanating from a sphere centered at the origin. Heyman 
|H89) translated Norris' result into the time domain using the analytic-signal 
(positive-frequency) Fourier transform. Norris and Hansen subsequently 
generalized the result to arbitrary bounded sources, both in the frequency 
domain |NH97 ] and time domain |HN97j . 



However, the representations |NH97| IHN97j are very different from (43). 
They express the weights of the complex-source beams in terms of the 
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spherical-harmonic expansion coefficients of g^^ , which requires only the field 
and not its normal derivative. On the other hand, since each of these coef- 
ficients involves an integration of the field over the entire sphere, it is not 
possible to express the weight of the complex-source beam emanating from 



an in terms of the incident field at that point, as in (43). Hence the ex- 



pansions in |NH97j and |HN97j are nonlocal, and consequently they do not 
have a straightforward physical interpretation like the one above. 



An electromagnetic analog of (43) has been published in |TPB7] and used 
in [TPB7Aj to accelerate the method of moments. 



The representations (34) and (43) can be further generalized to surfaces 



S other than spheres. It need not even be assumed that the source disks 
represented by the points of the analytically continued surface S must be 
tangent to S. However, this more general analytic continuation is more 
difficult than extending a single parameter like R. It does not work for all 
'regular' surface^for which a real Huygens representation holds because the 
integral expression is not necessarily analytic in a sufficiently large domain. 
To obtain an analytic continuation for a surface S, it is necessary to ensure 
that (a) the integration avoids all branch cuts, and (b) the area measure of 
S, which involves a Jacobian, can be continued analytically. These topics 
will be considered in future work. 



4 Gaussian pulsed beams as Huygens wavelets 



Care must be taken when transforming (38) to the time domain because the 

(45) 



integrand can grow exponentially in lo. Letting 

( = Cr + Ce=£,-iV, C = Cr + Ce, ^1 = Vr - Ve , 



the exponential in ( 38 ) is 



(46) 



Setting a' = a (as it will be after applying da'a), ( 19 ) shows that the bounds 
on r], as h varies with tCg fixed, are 



-2a < < (1 — 7)a where j = \ 1 



(47) 



Tegular surface is defined by Kellogg |K67| : see also [HY991 Chapter 2]. 
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The upper bound of ij is therefore positive whenever ^ 0. This divides 
the sphere Sji into the subsets 

S^{Xe, Xr) = {Rfi : r]r > r]e} (48) 
S^{Xe, Xr) = {Rh : rjr < rje] 

shown in Figure |5j As indicated, these sets depend on Xe and Xr- We call 



X' 




■y 



Figure 5: Decomposition (48l with R = 5, a = l,Xe = (0, 0, 2.5), cc^ = (20,0,0). 



The dark region (closest to x^) is S'^ and the light region is S'^. 



the frontal zone and the rear zone of Sr for the given emission and 
reception points x^^Xr. Note that the maximal value ?7e = o is attained 
when a^e is in the direction of — n, i.e., 



Ve 



I -^e 



(49) 



which gives the weakest contribution. This is a result of the opposite orien- 
tations of the reception and emission disks. 



To obtain the time-domain version of (38) choose a signal g{t), multiply 
both sides by g{io), and take the inverse Fourier transform. Formally, this 
gives 

g{t - r) 



87r2 



dn 

CrCe 



do; e-^^(*-^)5(u;) 



(50) 
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where we have exchanged the order of integration on the right side, which is 
justified if the double integral converges absolutely. If g is real, it suffices to 
compute its positive-frequency component and then take the real part. The 
positive-frequency component of g{t) is called its analytic signal: 

1 r°° 

m = ^J^ d^e-^-*5(^), K^) = H{u;)giu) (51) 

where H(uj) is the Heaviside step function. Taking the complex conjugate 
and using the reality condition g{uj)* = g{—u;) gives the negative-frequency 
component. 



1 f°° 1 







m* = - I du; e-*5(-^) = ^ / e--*5(^), 
^TT Jo Zir J_oo 

hence 

g{t) = 2Reg{t)- (52) 
If g{uj) decays sufficiently rapidly as w ^ oo, then the integral 

1 r"^ 

~9{r) = i^ duje-'^^g{uj), t = t + is (53) 
2vr Jo 

defines an analytic function of the complex time r. The domain of analyticity 
depends on the decay properties of g and the value of s. Formally, the 
positive-frequency part of (|50|) is therefore 



git — r) o? ^ f dh ^, , ^, 



provided the integral (53) defining g{t — Q converges absolutely for all h. 
Of special interest is the impulse 

g{t) = 6{t) g{uj) = 1 => g{T) = s < 0. 

ZTTIT 

The integral converges to the Cauchy kernel for s < and diverges for s > 0. 
The choice g{t) = 5{t) is very attractive since 

^-^^^ = P{x,t) (55) 
r 

is the retarded wave propagator, the unique causal fundamental solution of 
the wave equation: 

□P(x) = {d"^ - V^)P(x) = 4tt6{x), X = {x, t). (56) 
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P represents the wave radiated by the point source 6{x) at the origin of 
spacetime. It is 'fundamental' because it generates the field radiated by a 
general source g through 

F{Xr)= j d'^XeP{Xr-Xe)g{Xe) BF {x) = 4tT q{x) . (57) 

Thus, if we could obtain a pulsed-beam expansion for P, this would im- 
mediately give a similar expansion for all radiation fields F. However, it 



turns out that the divergence of (53) for s > makes this task very difficult. 



Equation (54) requires 

~g{t-C) = 9{t-^ + iv) 

both when r]r < rje and rjr > r]e- Numerical experiments have shown that 



while (54) 'almost' works with the Cauchy kernel, there is always a small 
but critical failure interval T = \ti,tt\ where it fails to converge. 

Note that disks are ideal for radiating beams (hence we have dish antennas), 
and recall that each point on Sr represents a tangent disk of radius a. Thus 
it is reasonable to try constructing a compressed representation of radiation 
fields by boosting contributions from the frontal zone S^, where r]r > rj, 



e-i 



and suppressing contributions from the rear zone Sj^, where r/r < r/g. The 



main contributions to (54 ) then come from the frontal zone, and this justifies 
the name 'compression.' 

However, the Cauchy kernel does this too well: it not only boosts contri- 
butions from the frontal zone; it makes them infinite, thus destroying our 
representation. We shall solve this problem with an elegant regularization 
which behaves naturally with respect to spacetime convolutions. This is very 
important because Huygens' principle is based on spacetime convolutions, 
as we shall see. Let 

gait) = ^—--, d>0. (58) 

\/TTd 



This is the Gaussian distribution with standard deviation a = (i/\/2- Al- 
though it seems that generality is lost by specializing to gd, this is actually 
not the case because 

d^O ^ gd{t)^5{t). (59) 

Therefore every continuous signal can be expressed as the limit of a super- 
position of translated versions of gd'. 

/oo 
dt'g{t')gd{t-t'). 
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That is, gdit) and its translates form a generalized 'basis' for signals. Define 
the Gaussian wave propagator 

P,(^) = Mi^, x = {x,t), r = \x\. (60) 



By (59), Pd converges to the retarded wave propagator as d ^ 0: 



lim Pd{x) = = p(^). (61) 

d^o r 

Its source is a 'Gaussianized' version of 6{x): 

aPdix) = 47rgd{t)5{x) = 47r5dix), lim 5dix) = 6{x). (62) 

d^O 



Just as P generates all radiation fields -F by (57), so does Pd generate their 
Gaussianized versions: 

Fd{Xr) = d'^XePdiXr - Xe)Q{Xe), lim Fd{x) = F (x) , (63) 

whose source is a Gaussianized version gd of g: 

DFd{x) =4tt d^Xe 6d{xr - Xe)g{xe) = 4:TTgd{x), lim gd{x) = g{x). 

J 

The Fourier transform of gd{t) is 

9d{u^) = e~''-y\ 

thus 



1 



oo 



9d{t) = ^l da;e-^-*e-'^'-'/^ (64) 

Both sides extend analytically to the whole complex time plane, giving a 
Fourier representation of the entire-analytic function gdi^)'- 



fl'd(T). 



vr d 
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The positive-frequency part of gd is 

27r Jo 27r Jq 



Thus 



Ur) = \ erfc (ir/d)g,(r) = "^^^^ (66) 



where 

erfc (ir/d) = ^ / due"" = 1 - erf (ir/d) 

is the complementary error function and w is the weh- known Faddeeva func- 
tion. Since both g^, and erfc are entire, so is ga- Define the function 

Hd{s - it) = Hd{-iT) = ^ erfc {ir/d), 

so that 



gd{T) = Hd{-iT)gdiT). (67) 



As ihustrated in Figure[6| erf (s/d) is a smoothed version of Sgn(s): 

erf (s/d) = / due " ~ Sgn(s), hm erf (s/d) = Sgn (s) 

and the smoothing is of order d, meaning that 

s < -d ^ erf (s/d) PS -1 and s>d ^ erf(s/d)f«l. 

Since 

1+ Sgn(g) ^ ^ /l, s>0 
2 \o, s<0 

is the Heaviside step function and 

2Hd{s-it) = 1 - erf (ir/d) = 1 + erf ((s - it)/d), 

Hd{s — it) is the analytic continuation of a smoothed version of H{s) with 

\imHd{s)^H{s). 
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Figure 6: The error function erf {s/d), here plotted with d = 2, is a smoothed 
version of Sgns to order d. As d 0, erf (s/d) — > Sgn (s). 



Again the smoothing is of order d: 

s < -d ^ Hd{s) ~ and s > d ^ Hd{s) « 1. 



(68) 



For small d, |i^d(s — it)\ is remarkably close to H(s) when \s\ > \t\. This 
can be seen in Figure [7j 

Furthermore, since 

erf (ir/d) + erf {-ir/d) = ^ erfc (ir/d) + erfc {-ir/d) = 2, 
Hd extends analytically the partitioning property H(s) + H{—s) = 1: 

Hdi-iT) + Hd{iT) = l. (69) 



Since ^^(t) is even, ( p7| ) and (69) imply 

gd{T) +gd{-T) = Qdij). 

But 



gd{i + is) 



vr d 



-2ist/(f 



9dit)e 



(70) 



(71) 



hence gd{t + is) grows exponentially when |s| > |t| and decays exponentially 



when |t| > \s\. The factor Hd{—iT) in (67) suppresses the negative cone 
s < —\t\, thus making gd{t + is) small everywhere outside the positive cone 
s > \t\. This is borne out in Figure pj 
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Figure 7: Plot of \Hd{s — it)\, shown from the side (left) and from below (right), 
where it is seen to be an approximation to a smoothed version of H(s) for \s\ > \t\ 
and have exponential growth for \t\ > |s|. The smoothing is of order d and the spikes 
along s = |i| are zeros. The phase of Hd is color-coded on the surface representing 
its modus IP9I. 



More precisely, the continuous-fraction expression |AS70| 7.1.4] for erfc im- 
plies that 



|t| ^ oo, s < =^> gdir) ~ 



1 



hence by ( 70 1 



|t| oo, s < 
The substitution r — r gives 



1 



2TriT 



|r| ^ oo, s > ^ Qdir) ~ gair) + 



1 



27rir ' 



and the two estimates can be combined into one that will be very useful^ 

1^ 



|r| ^ oo Qdir) ~ H{s)gd{T) + 



2TTiT 



(72) 



The 'small' value of g(i(r) in the region s < \t\ for large |r| is therefore the 
Cauchy kernel. Note that 



9d{-T) 



1 f°° 1 



2tt 



27r 



duo e 



-iiAiT-cp-up-lA 



(73) 



Equation (72 I is valid for s = since gd(i) ^ as |t| ^ oo. 
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Figure 8: Plots of \gd{t + is)\ (left) and \gd{t + is)\ (right) with d = 1. gd grows 
exponentially in the double cone \s\ > \t\ and decays exponentially in the double 
cone \t\ > \s\ , while gd grows exponentially in the single cone s > |t| and decays 
elsewhere. The dimples in gd are zeros of erfc(iT/d) (see Figure [?]) and the phases 
of gd and gd are color-coded on the surfaces representing their moduli |P9| . This 
shows the oscillation at the compression frequency p4| ) in the plot of gd{T) and its 
perturbed version (due to the complex factor Hd{—iT)) in the plot of gd{T). 



is the analytic continuation of the negative-frequency part gd{t)* of gd{t), as 



is also clear from (70) 



Equation (67) is remarkable. It shows that Hd{—iT) projects out exactly 
the positive-frequency part of ^^(t) by multiplication in the complex time 



domain, precisely as does H{uj) through (51) in the frequency domain. 



2x10' 
IxlO' 


- 1 X 10' 
-2x10' 




2x10' 
IxlO' 


-IX id' 
-2x10' 




Figure 9: The real part (left) and imaginary part (right) of gd{t + is) with d = 2 
and s — 8, together with their envelopes. 



Figure [9] shows the real and imaginary parts of gd{t + is) as functions of t 
with given values of d and s. They are very similar to those of the real and 



imaginary parts of gd{t + is) (71 ). 
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With g = ga, the positive- frequency analytic Huygens relation (54) converges 
absolutely: 



gd{t -r) q2 f dn 



(74) 



By (52), the Gaussian wave propagator (60) is given by 

„2 



D / X 9d{t-r) „^ /«% /" dn 
P.(x) = ^— = 2Re|-a.,.y — 



Ut - C) • 



(75) 



Carrying out the differentiations in (74) gives an expression more suitable 
for computations: 



gd{t - r) 



dn 



4^7 CrCe 



Se Vr 

Ce Cr 



+ (Ce - C)dt 



Ut - C) (76) 



where Cr = daCr and Ce = c^aCei as in (40) after setting a' = a. The 

:e 

and drgdir) 



derivative drgd{t — C) is easily computed. Since 



1 f°° 2 6' 

drHdi-ir) = —dr / dne-" = 



2r 

'd2 



5d(T"), 



we have 



drgdir) = gd{T)drHd{-iT) + Hd{-iT)drgd{T) 



2t 



9d{r) 



1 



(77) 



Note that 



Irl — > CO with Itl > s 



drUr) = 0{t-^) 



because the Cauchy kernel is canceled by ( 72 ) and the next term in the 



asymptotic expansion of ^^(t) is 0{t '^). Inserting this into (76) and us- 



ing (40) gives an expression without any derivatives, ideal for numerical 



computations. 



We shall now interpret (75) as a representation of Pd by a sum of pulsed- 



beam wavelets radiated by the disks I?(an) tangent to the sphere Sn. It 
suffices to work with the positive- frequency part ([74|). Recall that 



•JUt' JU fi 
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represents the spacetime 4- vector from the emission event Xe to the reception 
event x^. Consider the intermediate complex event given bjj^ 

z = (a'n, r), where r = te + Ce (78) 

is the emission time te plus the complex travel time C,e from x^. to afi. Define 
the Gaussian pulsed-beam propagator from z to Xr by 



Fd[Xr - Z) = = , t = tr - te, C = Cr + Ce- 



(79) 

This represents the complex wave amplitude radiated by a'h at the complex 



time r and received at at time tr- Thus (74) reads 

,2 



f ~ 1 

Pd{Xr - Xe) = — da' a / dh Pd{Xr - z) — . (80) 

4vr J Ce 



The general Gaussianized solution Fd{xr) in (63) is therefore given by 

Fdixr) = 2Re F(i{xr) where 
~ Oi^ (f~ o( X ) 

Fd{Xr) = da' a J '^'^ J dXePd{Xr " ^)~^- (^^) 

More will be said about pulsed-beam representations of general solutions in 
Section [7l 

We now show that Pd{xr — 2:) is a pulsed beam radiated by T>(a'h) which 
propagates along n, with the propagation along —h suppressed by gd as in 
Figure [s] The factor C,~^ represents the attenuation suffered by the spherical 
wave emitted by the point source at Xe while propagating to ah. Thus we 
have a picture, shown in Figure [TO] of a spherical wave emitted at Xe and 
received at ah with 'reception amplitude' C~^, then re- radiated from a'h 



as a pulsed beam and finally received at Xr ^^ The idea of analytically 
continued fields as reception amplitudes by complex-source disks will be 
explained in more detail in Section |8| 



By (79), 



P.(x.-.) = MizO= ^^(^-/ + -^) . (82) 



^Recall that ("e = and i^,- = C'-(q') and we set a' = a after applying da'a 

dec' — da- 

^°We are ignoring the derivatives da'a, so this interpretation is somewhat schematic. 
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Figure 10: The factor in (80 1 represents the reception ampHtude at ah due to 
the attenuation suffered in propagating from to an, and Pd{xr ~ z) represents 
the propagation of a pulsed beam from the disk T) to Xr- 



The properties established for gd show that the magnitude of Pd[xr — z) is 
an increasing function of r/ that attains its maximum values in the frontal 
zone nerest to Figure 5l Some insight can be gained by noting 

that 

and expanding 

Pd{xr -z) = C^'Hdirj - i{t - e))e'''/^'<7d(t - i)e-^'^^'-^'^'''\ (83) 

• At a given time t, the factor gfi{t — ^) ensures that Pd{xr — z) is 
concentrated on a shell of thickness ~ 2d around the surface = t. 
Pd has significant values only when t is in the range of ^, which varies 
with h £ S'^ over a positive interval containing the line of sight time 
r = \xr — Xe\, the minimum time required to travel from to Xr at 
speed c = 1. If instead we vary Xr = {xr,tr) but fix Xe and h, this 
means that the oblate spheroid given by 

is a wavefront of Pd{xr — z) expanding with t = tr — te- This gives 
a direct meaning to ^r- it is a variable whose level surfaces are wave- 
fronts. 
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From the behavior gd, it follows that the factor Hii{r} — i{t — £,))e 



in (83) is an increasing function of t] that boosts the incoming wave 



when 7/ > while suppressing it when 77 < 0. 

• Due to the factor e~'^*^^*~^^/'^^, Pd{xr — z) oscillates at the compression 
frequency 

Mv) = ^ , (84) 

which depends on Xr for given h and on n for given x^. This is 
perturbed slightly by the phase of Hd{ri — i{t — ^)). 



By interpreting every factor in (83), we have thus understood Pd{xr — z) as 
a pulsed beam with wavefronts O^^ propagating along the semi-hyperboloid 
TLri^ at the compression frequency cod- 



Consider the limit of (80) as a, a' 0: 

I' ~ 1 

Pd{Xr -Xe) = — Or'r / dh Pd{Xr -x) — , 

where x = {R'h, te + rg) is a reception event on Sr' at the arrival time te + rg 
of a spherical wave radiated from x^ at te- As a function of x,., Pd{xr — x) is 
the positive- frequency part of a real Huygens wavelet emitted from By 
complexifying the sphere, we have deformed the original spherical Huygens 



wavelets to pulsed beams ( 79 ) : 

Pd{Xr -x) ^ Pd{Xr - z). 

This deformation acts on space so that spheres become oblate spheroids 



and cones become semi-hyperboloids, as in (16). In the process of being 



deformed, the spherical Huygens wavelets are compressed in the forward 



direction and stretched in the backward direction. Being complex, the 



Although Pd{xr — x) = 2 Re Pd{xr — x) is not a wave because gd{t — r) does not 
oscillate, applying the derivative da'a gives it some oscillation. For example, 

dc^gdit - C) = -CdtgS - C) = - C)5d(t - C) 
is a one-cycle wave. 

^■^In a certain sense, they are Doppler scaled positively in the forward direction and 
negatively in the backward direction [K94| . 
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compression introduces a phase which gives a measure of its strength. This 
is why we cah uod the 'compression frequency.' Note that 



as expected. 



The time-domain radiation pattern of a radiation field F{x, t) with cyhn- 
drical symmetry is, by definition |HY99| . the function T{0,t) satisfying the 
far-field relation 



F{x,t) 



J^{e,t-r) 



To compute the radiation pattern of Pdi^r — z) relative to the coordinate 
system of the disk T>(a'n) , assume the observer is far from the disk, ^'^ Taking 
a' = a for simplicity, (23) gives 

Vr ^ a ^ ~ l^r, rjr ~ aVr ' ^ = a COS Or, 

so that 

^ ~ + ^e, r/ ~ a cos Or - Ve- 



The factor (r^ in ([79]) can be approximated by ^ since Vr ^ a and < 
re < 2a. Thus 

Fd[Xr — ~ where t = t — t,e- 

rr 

Hence the radiation pattern of Pd{xr — z) is 

T{9r,t) = ~gd{t' + i{acos9r - r^e)). (85) 

The peak radiation time is t' = 0, when t = is the arrival time at 2?(an) 
of the emitted wave. Figure 11 shows polar plots of the peak-time radiation 
patterns for two values of a and the the two extreme cases with 



r/e = 7oa and rje = a. 



where 



7o = luTi 7 



i?2 ^ 

1 n ^ = , < 1 

+ a? ^/WTa? 
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i?=10 a = 5 d = 2 
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Figure 11: Peak-time radiation patterns of Pd{xr — z) for d = 1,2 and a = 5,50. 
In each case we have plotted the beams with the weakest pattern {rje = a) and the 
strongest pattern (rje — 7Qa). For a — 5 and d = I, the weakest pattern is so weak 
that it cannot be seen. 



as in ( pO| ). This lower bound applies to every source supported in Sr. 

Since \gd{t + is)\ is an increasing function of s, the upper bound ije = a is 
expected to produce a weaker pattern than the lower bound rj^ = 70a, as 



already discussed beneath (49). This is borne out in Figure 11 where the 
pattern with r]e = a for a = 5 and d = 1 is so weak that it is invisible. On 
the other hand for a = 50 the disk is so large as to dwarf the sphere. Since 
7o ~ 0.98 in this case, there is not a great difference between the lower and 
upper bounds of r/g. This is the reason why both patterns are visible in the 
lower figures, and why both are much weaker than the patterns with a = 5 
and r]e = 7oa- 

We have thus established (|80[ ) as a pulsed-beam representation of Pd{xr—Xe)- 
Since the pulsed beams Pd{xr — z) are deformations of Huygens' spherical 
wavelets, it is reasonable to call them pulsed-beam Huygens wavelets. By 



taking the real part and convolving with a general source g, as in (63), we 
obtain a pulsed-beam representation of the Gaussianized version Fii{x) of 
general solution F{x). 



Since we are keeping n fixed but varying Xr, it is unnecessary to assume that S> R- 
Only the relative vector — x,- — Rn enters the above discussion. 
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5 Some remarks 



1. The radiation patterns of extended sources generally have sidelobes, 
which are interference patterns between parts of the wave arriving 
from different parts of the source. Sidelobes of beams often stray 
widely from the intended direction of propagation, causing problems in 
applications such as communications, remote sensing, and radar 



However, note that the radiation pattern ( 85 ) is real at the peak time 



t' = and it decays monotonically with increasing Or, as confirmed by 



Figure 11 It therefore has no sidelobes, and that makes it potentially 
very useful. If we include the time-dependence around t' = 0, the 
radiation pattern acquires a phase factor e"*'^''*' and 2 Re J='{9r,t), 
the radiation pattern of Pd{xr — z), acquires sidelobes. But these are 
confined to the narrow envelope of 2\J^(9r,t)\ and do not cause the 
usual problems. 

2. To fully justify the name 'pulsed-beam propagator,' consider Pd^Xr — z) 
as a function of Xr with z fixed. It is singular on the branch cut 'D(a'h) 
of Cr(a') and analytic elsewhere, hence 

Xr V{oln) UrPdiXr - z) = 

where Or is the wave operator with respect to Xr- Pd{xr—z) is therefore 
the wave radiated by the disk 'D{a'h) . The precise source of Pd{xr — z) 
is a generalized function 5d{xr — z) supported on Xr G P(a'n); see 



Equation (108) in Section 8 



3. Taking the complex conjugate of (74) and substituting a, a' — > a*, a'* 



(which is permitted since the left side is independent of a and a' and 



the domain (33) is symmetric under conjugation) gives the negative- 



frequency component in the form 

~gd{t - r)* f dh 



, / ~g^(t-C)\ (86) 

47r J CrCe 



where we have used the reality conditions ( |18[ ) 

Ce(a*)* = Ce(«), Cr(a'*)* = Cr(a'). 

Thus gd{T) satisfies the reality condition 

9d{T*Y = gdir), 
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which also foUows directly from ( 65 ) . Adding (1741 and ( 86 ) gives an 



alternative form of the analytic Huygens relatior 



9d{t 



a 
47r 



5„ 



dn 

CrCe 



9d{t-0, 



(87) 



which is simpler than ( 75 ) as it does not split up the positive and 



negative frequencies. However, we find that while (87) is numerically 



valid, it does not lead to a compressed representation of radiation 



fields. The problem is the substitutions a 
R + ia with a > 0, 



a , a 



a 



For a 



Rh 



lan. 



Hence the disk P(a*), while still tangent to Rh, radiates a pulsed 
beam along — n, i.e., to the interior of the sphere. Eventually, this 
beam leaves the sphere and continues to propagate, weakened, in the 
direction of — n; but this is clearly an inefficient way to represent 



radiation. Although (87) is mathematically correct in the sense that 
the integral on the right converges absolutely to Pd{xr — Xe), this 
inefficiency shows up in the appearance of very large numbers which 
spoil the compression and easily overwhelm computational software, 
thus introducing huge errors; see the discussion at the end Section [9j 

4. In view of the previous remark, we can say that the positive-frequency 
part of (50) is 'good' while its negative- frequency part is 'bad.' The 



situation would be reversed for the interior problem, where a source 
is given outside of Sr and we seek to represent the field inside Sr 
as a superposition of pulsed beams. The pulsed-beam analysis and 
synthesis of interior fields is very similar to that of exterior fields and 
will be treated elsewhere. 



5. The pulsed-beam representation (81) of general radiation fields sug- 



gests an important application: given a receiver at Xr, the most signif- 
icant contributions are expected to come from disks radiating approxi- 
mately in the direction of Xr, whose centers Rh are in the frontal zone 
S^. That is, by using only the 'relevant' wavelets propagating toward 
a given observer, we obtain a compressed representation of Fd{xr). 



This is discussed in greater detail in Sections^ and 10 



Equation (87) can also be obtained directly from (50 1 with Qdi^^) ~ e 



V4 
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6 Huygens reproducing relation for pulsed beams 



The time-domain version ( 80 ) of the analytic Huygens principle treats emis- 
sion and reception asymmetrically: propagation from Xe to z is represented 
by Ce^: whereas propagation from z to Xr is represented by Pd{xr — z). In 
this section we construct a more complete picture of this process which has 
a detailed and appealing physical interpretation. For this we shall need to 
Gaussianize both the emission time te and the reception time tr- Thus let 
de and dr be Gaussian duration parameters for te and U and let 

d=^/di+dl 

which will be the duration parameter for the entire transmission process. 
Let 

Za = {an,T), Zai = (a'n, r), where t = t + is (89) 

is a free complex time variable. When r = tg + Cej t = te + £,e and 
s = r]e, (89) reduces to (78). The propagations from Xe to z^ and z^' to Xr 
are governed by 

5d,(T -te- Ce) 



Pde^Zd Xe) 
Pdr{Xr - Za') 



Ce 

gdXtr - T- Cr) 
Cr 



with 



1 

~9dAr -te-Ce) = ^ ^^^-^^{t-U-C.+^s)^-dWl^ 

1 

9dAtr -T-Cr) = :^ du; ^-i^{tr-t-Cr-is)^-dWl\ 

Applying the Fourier transform in = t — te to the first equation and in 
t" = tr — t to the second equation gives 

/oo 
dt' e''^''gdSt' -Ce + is) = H{uj) ^^-^(.(^^-i^) ^-^W 1^ (90) 
-00 

/oo 
dt!' e^^'"gdr{t" -Cr- is) = H{uj) e^^(f'-+^^)e-'^''^'/^ 
-00 
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Multiplying (pSj) by i7(a;)e-'^'^'/^ = H{ujfe-'^^'^'/^ gives 



H{uj 



Cr 



Ce 



a 
47r 



3iaj(Ce-is)g-c(^tj2/4 
Ce 



Taking the inverse Fourier transform and writing the time variable as tr — te 
gives 

gd{tr -te-r) 



dn'n I dh I dt 

-oo 



gdriir - T -Cr) QdS'^ -te- Ce) 



r An 
We have thus proved the following result. 



Cr- 



Ce 



Theorem 2 The Gaussian pulsed-heam propagator satisfies the following 
complex spacetime Huygens reproducing relation: 



a 



Pdi^r Xe) — , da'a 

47r 



dn 



dtPd^{Xr - Za')PdcXZa - Xe). 



(91) 



Remark 1. The complex spacetime 4-vector Za = {ah,t + is) represents 
a pulsed receiving disk with a 'Gaussian' reception interval [t — de,t + de], 
which we denote b}j^ 

Vd,{z^) = V{ah) X Id^t) C where Id^t) = [t - 4,t + 4]- (92) 

Just as ah G C'^ represents the extended object D(an) in space, so does 
z G represent the extended object 'Dd^{z) in spacetime. Similarly, Za' 
represents a pulsed emitting disk 



VdXza')=V{a'h) X hXt)- 



(93) 



The relation between pulsed-beam emitters and receivers will be explained 
in greater detail in Section Isl 



Remark 2. Equation has a simple interpretation which, unlike (80) 



treats emission and reception symmetrically. It states that the spherical 
wave emitted from the point source Xe is received by ^^^(zq), then im- 
mediately re-emitted by Vd^iza'), and finally received at Xr- The direct 



'The role of s = Im r will be explained below. 
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propagator Pd{xr — Xe) is recovered by integrating over all directions n and 
intermediate times t and then applying (a'^ /Air) da' a- 

Remark 3. The integral over t can be viewed as a contour integral in r, 
with the left side independent of s due to analyticity. In fact, s is allowed 
to depend on n (and even on t). For a general emission source ge{xe) and 
receiving source Qr{xr), where there are additional integrations over Xg and 
Xr, s may also be allowed to depend on Xg and Xr- 



Remark 4. The formal symmetry of (91) with respect to emission and 
reception is of more than purely academic interest. To formulate pulsed- 
beam representations for the interior field given an exterior source g, we 
must reverse the roles of Xe and Xr and use the advanced wave propagator: 



P{x) P'{x) 



6{t + r) 



(94) 



Then (91) transforms to a pulsed-beam representation of the interior field 
but (80) fails to do so. However, to get an efficient representation, we must 
also replace a, a' by their complex conjugates, as discussed beneath 
since we now want the pulsed beams to propagate inward. 

The properties of gd place some practical constraints on s. Let 



Tr 



T 



Ce 
Cr 



t - 



te- ie + «(s - 
-t - ir + iiVr 



(95) 



s), 



Pdei^Ol Xg) 



(96) 



so that 

^ and hX-r-Za,)=~^. 

The compression frequencies of the interior and exterior pulsed beams, de- 
fined as in (84), are 



2(g - rje 
<R 



UJr 



dl 



(97) 



and the propagators in (91 ) will be very small unless 

UJe > 0, UJr > 0, 

respectively. If both inequalities hold, they imply 

ije < s <rir 



(99) 
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which is consistent with r] = r]r 
r]r + Ve 



rje > 0. For example, 



and 



dr 



(100) 



Since rjr > rje for the dominant beams, the condition (98) is indeed satisfied 



by (100), showing that (99) is sufficient as well as necessary. Although our 



proof of (91) is theoretically valid for all choices of s, the computation can 



be expected to be inefRcient for values of s violating (98). For example, if 
s is large and positive, then Pd^{za — Xe) is very large and PdX^r — Za') 
is very small. Conversely, choosing s large and negative makes Pd^ very 
small and Pd,. very large. Such choices introduce unnecessary noise into the 
computation, thus reducing its efficiency and even causing errors when the 
machine capacity is exceeded, which does in fact occur rapidly due to the 



exponential growth. Note that choosing s in the 'good' interval (99) makes 
it dependent on n, which is permissible as explained above. 

Numerical calculations confirm that all values of s in the interval [r/e, rjr] give 
stable results and that instabilities build up rapidly when s strays outside 
this interval. 



7 Analytic Huygens relation for general solutions 

Let g{x) = q{x, t) be a time-dependent source distribution bounded in space, 
and choose R so that Q{x,t) is supported in the open ball \x\ < R at all 
times, The radiated field F is given by 



F{x) = j dXePix - Xe)g{Xe) 



(101) 

where P is the retarded wave propagator ( |55| ). Now Gaussianize the emission 
time te by 

/oo 
dt'^gd,{te - t'jQ{Xe,t'^), de > 0. 
-oo 

Then the wave arriving at x is (by the associativity of convolutions) 

P'deix) = j dXeP{x - Xe)QdSXe) = j dXe Pd^{x - Xe) QiXe) (102) 

^^This is always possible if g is compactly supported in time as well as space. If g 
is spatially bounded at all times but does not remain in a bounded spatial region (for 



example, if it drifts at some velocity v 0), then a generalization of (911 based on a 
spacetime version of Green's second identity must be used. 
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where Pd^ is the Gaussian propagator ( 60 ) . The spatial integral is over the 



support of Q in the interior of the sphere Sr. Both sides of (102) can be 
continued analytically to the complex spacetime points 

X ^ Za = [ah., t + is) 

since the integration over x^, does not encounter any of the branch cuts of 
C,e = \/ {otfi — Xf,)"^. The analytic continuation of the positive-frequency part 
of dTol is 



FdS^a) = I dXePdS^a- Xe)Q{Xe). 



(103) 



Convolving (91) with q gives 



a 



Fd{Xr) = ^da'a / dn / dt P^^Xr - Zai)FdSZa) 



(104) 



where the reception time tj. has also been Gaussianized by convolving with 
Qd^^ so that the total duration parameter \s d = \/d1 + d^.. 



The analytic Huygens relation (91) for propagators thus implies a pulsed- 



beam representation for arbitrary solutions with spatially bounded sources. 
As will be explained in Section|8j the coefficient Fd^{za) in this superposition 
is the reception amplitude of the interior field Fd^ by the pulsed disk VdSza)- 



8 Pulsed-beam reception and emission 



The purpose of this section is to justify the interpretation of Fd^ {za) in ( 104 ) 
as the reception amplitude of the field F by the pulsed disk Vd^iza) defined 



in ( 92 ) . By the wave equation 

□F(x) = 4TTg{x) 



(103) can be written as a relation between Fd^Za) and F{xe), 

1 



Fd^Za) = ^ J dXePdAZa - Xe)OeF{Xe), 

where Dg is the wave operator in and 

Za = {ah, t) = X + iy, x = {Rh, t), y = {ah, s) 



(105) 
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Integrating by parts twice give; 



T7 



1 

47r 



(106) 



To make sense of this, note that in real spacetime {y 0) we have ( |62[ ) 

□ePd,(rC - Xe) = ^T^Qd^it - te)S{x - Xe) = 4TT5d,{x - Xe), (107) 

whose positive-frequency part is 

OePdS^ - Xe) = 4:7rgd,{t - te)d{x - Xe) = 4:7rSdSx - Xe). (108) 

We now define the source distribution 6d^{za — Xe) of Pdei^a — Xe) by ex- 
tending this to complex spacetime: 



'^■K^dS^a - Xe) = UePdS^a - Xe). 



(109) 



Since PdS^a — Xe) is analytic whenever Xe ^ I?(an), it follows from ( |107 ) 
that 'dePdS^a — Xe) = at such Xe. Hence the distribution 5d^{za — Xe) is 
supported in Xe G T>{afi) at all times. It is also localized in a 'Gaussian' 
sense in time around the interval Id^it), thus it is effectively localized in the 
pulsed disk Xe G Vd^iza). 

In other words, the wave operator De ignores all the analytic behavior of 
Pdei^a — Xe) and nails down its singular behavior, consisting of a discon- 
tinuity across the branch cut P(an) and infinity along the branch circle 
&D{an). All this is possible only in the distributional sense. Combining 



(106) and (109) gives 



Fd^Za) = / dXeSd^Za - Xe)F{Xe) 



(110) 



which is the analytic deformation of F. Equation (110) shows that Fd^ {za) 
is the reception amplitude of F(xe) by the receiving source Sd^{za — Xe), 
confirming our claim. 

Remark. The distribution 6d^ includes a dipole layer, represented by a 
first-order differential operator acting on F{xe) |K0| IK3j . Consequently, the 



right side of (110) contains the values of both F and its partial derivatives. 



^"^ Pd^{za — Xe) is singular when cce G 'D{an), so the right side of ( |106[ l must be treated 
carefully. The wave operator De acts on Pd^Za — Xe) in a distributional sense, just as it 
acts on Pd^ {x — Xe) to give 'in5d^ {x — Xe). The resulting distribution Sd^ can be computed 
rigorously using the methods developed in [KOI IK3I IK4I IK5I ID8| and will be studied in 
detail elsewhere. Here we explain the main ideas in an intuitive and informal way. 
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9 Compressed representations of radiation fields 



We now demonstrate that the computational properties of the Huygens rep- 
resentation depend strongly on the disk radius a and there can be a signifi- 
cant advantage to choosing a complex sphere over a real sphere in numerical 
calculations. 

Throughout this section, the source is at = (0,0,2.5) and has Gaussian 
time dependence gd{t) with d = 0.3\/2- The reception point is in the far 
zone on the positive x-axis. We consider three spheres with R = 10 and 
disk radii a = 0,5,50. For a = 0, the sphere is real and the pulsed-beam 
representation reduces to the classical Huygens representation. 

Let 



2 Re 



a 



47rCrCe 



+ (C - C)dt 



~9d{t-0}, C 



(111) 



where the dependence on Xe and Xr is implicit and dtgd{t — Q is given by 



(771. Then (76) reads 



Pd{x,t) = 2Rq Pd{x,t) = dhI2{h,t), {x,t 



X 



(112) 



Figure 12 shows evaluated at the peak time t = r, where Pdi^, t) attains 
its maximum value gd{0)/r. The top plot shows that for the real sphere 
is significantly nonzero only on a ring centered on the line of sight point 
from Xe to Xr- We can explain this behavior as follows: for a = 0, 

Ce = {Rn — Xe\ = Ve and C,r = |a^r — Rn\ = Tr 
satisfy the triangle inequality 

y I ^6 f j tl/^ tC g I J 

with equality if and only if Rfi is the line-of-sight point. Now 

~9d{'r - C) = Hd{-i{r - Vr - re))gd{r - Vr - re) 

decays as a Gaussian m r — — perturbed by Hd{—i{r — r^ — re)). Com- 
bined with the effect of the derivatives coming from da'a^ this perturbation 
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Figure 12: Magnitudes of with i? = 10 at the peak time t — r and a — Q (top), 
a = 5 (middle), and a = 50 (bottom). The emission point is x^. = (0,0,2.5) and 
the reception point is in the far zone {\xr\ S> |a|) on the positive x-axis, outside 
the figure. The line of sight from Xe to Xr is indicated by the long arrow. 
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Perspective view 



Side view 



z 



Spherical cap 




Figure 13: A spherical cap with cap angle (i centered on the x-axis. 



displaces the maximum from the line of sight (r^ + re = r) to a small circle 
centered at the line-of-sight point. The derivatives cause to vanish at 
the line-of-sight point and oscillate near the circle, thus creating the ring 
pattern seen in Figure [T2| 

For the two complex spheres in the middle and bottom plots, g^it — 
suppresses the points Rfi with t/.^ < r](. and boosts those with r]^ > rje in the 



frontal zone S]^ of (48). Simultaneously gd{r — S,) suppresses points with 



|r — ^1 > d, and the winners of this tug of war are the points in the small 



spot centered at Rxr in the middle and bottom plots of Figure 12 This 
spot becomes more and more concentrated near Rxr with increasing a, as 
shown in Section [Tol 



It is important to note that for large a is concentrated around the point 
Rxr nearest to the receiver regardless of the location of the emission point. 
This means that 22 large a would remain concentrated near Rxr, as 



in Figure 12 even if we replaced the point source at Xr with an arbitrary 
volume source g{xe) supported throughout the interior of Sr. This can 
be understood by noting that a disk source becomes more directive with 
increasing radius, hence fewer disks are required to achieve a given accuracy. 

To illustrate the advantages of using complex spheres with large value of a, 
we investigate the accuracy of the field at Xr obtained by including only the 
pulsed beams radiated from a spherical cap centered around Rxr with maxi- 
mum angle f3, as in Figure [T3l For the parameter values under consideration, 
/5o = 15° represents a cap where the line of sight from x^ to Xr just grazes 
the upper edge of the cap. We compute the maximum error over all time 
of the field at Xr, relative to the maximum value gd{^)/T~- This calculation 
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Figure 14: Compression: The maximum error at a far-zone reception point on 
the positive x-axis as a function of the cap angle /3. The emission point is on the 
z-axis at z = 2.5. The three spheres have R — 10 with disk radii a = 0, a = 5, and 
a = 50. 



is designed to simulate the realistic situation where the time dependence of 
the source is unknown and the source may emit a series of pulses that are 
spread out over time. In particular, the computed error bounds also hold 



for time-harmonic fields. Figure 14 shows this maximum error as a function 
of /3 > /3o for the three values of a. A dramatic error reduction is obtained 
by increasing a. For example, the errors e{a) for P = 45° are 

e(0) = 27.9%, e(5) = 8.3%, e(50) = 1.5%, 

a reduction of nearly 20:1 when using the pulsed-beam representation with 
o = 50 compared to the real Huygens representation! This has important 
practical implications in numerical calculations, where a certain error level 
must be often achieved. For example, assume that the error of the field at 
the reception point must be less than 2%. Then the required cap angles 
/3(a) are 

/3(0) = 152°, /3(5) = 89°, /3(50) = 38°. 

Hence, with a = 50 we need to include pulsed-beam wavelets over just 11% 
of the sphere, whereas a = would require spherical wavelets over 94% of 
the sphere. 
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Remark. Consider the efficiency of the 'alternate' expression using (87), 
whose negative-frequency component was obtained by letting a — > a* and 
cJ — > OL* in the conjugate pulsed-beam expansion (76) of Prf(a;,t)*. This 
gives the expression 

„2 



a' 



dh 

CrCe 



+ {C-Odt 



9d{t-C) ^ J dnj^{h,t), 

(113) 



where gd has been replaced by gd in the integrand. As noted under (87), 
this means that the negative-frequency pulsed beams propagate along — n 
instead of fi, thus traversing the sphere and spoiling the efficiency of (75). 
To demonstrate the enormous difference between (112) and (113), con- 
sider again the complex sphere with a = 10 -|- 50i, = (0,0,2.5), and 
Xr = (200,0,0) in the far zone. Since (113) does not have a factor like Hd, 
there will be significant contributions from the back side of the sphere. In 
particular, consider 



n 



-Xr 



(-1,0,0) C = Cr + Ce = 220.01 99.94i 
Choosing d = .3-v/2 as before and t = ^ = 220.01, we find 

9dit - C) = gd{-99Mi) = 1.336^^^89 ~ 10™. 



The contribution from the back point —Xr in (113) is seen astronomical! By 
comparison, the counterpart of gd{— 99. 94i) in (112) is 

gd{-99Mi) = 0.00159. 

This illustrates that (112) and ( |113 ) are on opposite sides of the computa- 
tional efficiency spectrum. 



10 Pulsed-beam representations for large a 



Figure 14 suggests that the compression keeps improving as a increases, 
and it is natural to wonder if there is an 'optimal' value of a beyond which 
the benefits of further increase will diminish. We therefore analyze the 
asymptotics of the pulsed-beam representation for large a. 

Define 6, 9e, and ip by 

cos 9 = Xr ■ n, cos 6e = ff, ■ h, cos 1p = Xr ■ Xe- (H^) 
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Suppose that R and that Xr is in the far zone (24): 

R <^ a and ja;,.! » |a|. 
Since \xe\ < R <^ \xr\, we have 



(115) 



r = \Xr — Xe\ = \Xr\ — Xr ■ Xf, + 0{R/\Xr\) ~ I^E^I — \Xe \ COS l/j- (116) 



Also 



and 



thus 



Ce = V {re + iafiY = re cos Be + 



la 



2a? 



+ OiR^a^' 



Cr = \J {xr — afiY = \xr\ — RcosO — iacosO + 0{a^ /\xj.\)^ 



C ~ Te COS 9e + \Xr\ — i? COS 9 + ia 



r1 s\v? 9e 
2a? 



This shows that for large a, the frontal zone ?] > is given by 



S^-^ iRh : sm{9/2) < 



Tp sm 



e oi±i ue 



2a 



(117) 



(118) 



As a increases, SX shrinks to the point {ii^^} nearest to the receiver. At 



the center of S'^ we have n = x^ and 



hence 



Tg = {RXr — Xe)"^ = R^ + |iCe|^ — 2i?|a;e| COS "0 
Te COS 9e = {Rx^ — Xe) ■ X^ = R — \Xe\ COS ijj, 



h = Xr =^ Tg sin^ 0e = |a^e|^ sin^ ^ = 6^, 



where b is the length of the projection of Xe to the plane orthogonal to Xr- 
Therefore 



n = Xr ^ T] 



2a' 



The real part of (117) gives 



h ^ Xr =^ ^ {RXr — Xe) ■ Xr + \Xr\ — R = \Xr\ — \Xe\ cos ^, 
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thus by ((TTel) 



n 



(119) 



We have thus established that C = i — if] behaves as follows as a increases: 



^ — > COS + I a^r I — R COS 9 and 
^ ^ r and r/ ~ /2a. 



-oo 



The asymptotic formula (72 ) shows that as a increases, gdit — C) 
a shrinking cap centered at Rxr while in the center we have 



outside 



n 



~gd{t-C)^~9d{t-r + ih^/2a). 



Figure [s] shows that the phase of gd{t + is) remains approximately constant 
for s < 0. Hence gd{t — C) decays smoothly without oscillations as i% moves 
away from Xr- Furthermore, it remains bounded throughout the sphere as 
a grows, and its maximum magnitude approaches \gd{t — ^)|- 

Thus the compression keeps improving as a increases, and there is no limit to 
how small the cap angle (3 can be for a given error limit. But as a increases. 



the sampling rate used in the computation of the integral (112) must also 
be increased to capture the rapid variation of gd{t — C) near n = x^. We 
shall investigate the numerical consequences of using large values of a more 
fully in future work. 
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